home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_3 / a3_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  4.3 KB  |  165 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 3.3 (PA = LU Factorization with Pivoting).
  9. % Section    3.6,    Triangular Factorizaton, Page 175
  10. echo on; clc; format long; clear
  11.  
  12. % This program solves linear systems AX = B.
  13.  
  14. % The method is LU factorization with pivoting.
  15.  
  16. % 1. Factor  PA = LU 
  17.  
  18. % 2. Get    LUX = PB
  19.  
  20. % 3. Solve   LY = PB  for  Y.
  21.  
  22. % 4. Solve   UX = Y   for  X.
  23.  
  24. % Remark. lufact.m and lusolv.m are used for Algorithm 3.3
  25.  
  26. pause % Press any key to continue.
  27.  
  28. clc;
  29. % Solve the system AX = B  where:
  30.  
  31. A = [1    5    4   -3;
  32.      4    8    4    0;
  33.      1    3    0   -2
  34.      1    4    7    2];
  35.  
  36. B = [-4;  8;  -4;  10];     % Enter B as a column vector.
  37.  
  38. [LU,row,det1] = lufact(A);
  39.  
  40. [X,Y] = lusolv(LU,B,row);
  41.  
  42. R = B - A*X;                % Verify that the residual is small.
  43.  
  44. pause % Press any key to continue.
  45.  
  46. [n n] = size(A);
  47. I = eye(n);
  48. for k = 1:n,
  49.   P(k,:) = I(row(k),:);
  50. end
  51. L = I;
  52. for c = 1:n-1,
  53.     L(c+1:n,c) = A(row(c+1:n),c);
  54. end
  55. U = zeros(n,n);
  56. for c = 1:n,
  57.     U(1:c,c) = A(row(1:c),c);
  58. end
  59.  
  60. Mx1 = 'Implementation of the LU factorization method.';
  61. Mx2 = 'The matrix is A =';
  62. Mx3 = 'The vector B is displayed as B` =';
  63. Mx4 = 'The solution to AX = B is displayed as X` =';
  64. Mx5 = 'The permutation matrix is P =';
  65. Mx6 = 'The lower-triangular matrix is L =';
  66. Mx7 = 'The transformed vector PB is displayed as (PB)`';
  67. Mx8 = 'The solution to LY = PB is displayed as Y` =';
  68. Mx9 = 'The upper-triangular matrix is U =';
  69. Mx10= 'The vector Y is displayed as Y` =';
  70. Mx11= 'The solution to UX = Y  is displayed as X` =';
  71. Mx12 = 'The residual R = B - A*X is displayed as R` = ';
  72. clc,echo off,diary output,...
  73. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(A),disp(Mx3),disp(B'),...
  74. disp(Mx4),disp(X'),diary off,...
  75. disp('Press any key to continue.'),pause,diary on,...
  76. clc,disp(Mx5),disp(P),disp(Mx6),disp(L),disp(Mx7),disp((P*B)'),...
  77. disp(Mx8),disp(Y'),diary off,...
  78. disp('Press any key to continue.'),pause,diary on,...
  79. clc,disp(''),disp(Mx9),disp(U),disp(Mx10),disp(Y'),...
  80. disp(Mx11),disp(X'),disp(Mx12),disp(R'),diary off, echo on
  81. pause % Press any key to perform LU factorization with pivoting.
  82.  
  83. clc;
  84. % Solve the system AX = B  where:
  85.  
  86. % A is the matrix from the previous example.
  87.  
  88. % Notice that we do not need to enter it.
  89.  
  90. % The LU factorization is used and not recomputed.
  91.  
  92. B = [-5;  4;  -3;  7];      % Enter B as a column vector.
  93.  
  94. [X,Y] = lusolv(LU,B,row);
  95.  
  96. R = B - A*X;                % Verify that the residual is small.
  97.  
  98. pause % Press any key to continue.
  99.  
  100. [n n] = size(A);
  101. I = eye(n);
  102. for k = 1:n,
  103.   P(k,:) = I(row(k),:);
  104. end
  105. L = I;
  106. for c = 1:n-1,
  107.     L(c+1:n,c) = A(row(c+1:n),c);
  108. end
  109. U = zeros(n,n);
  110. for c = 1:n,
  111.     U(1:c,c) = A(row(1:c),c);
  112. end
  113.  
  114. clc,echo off,diary output,...
  115. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(A),disp(Mx3),disp(B'),...
  116. disp(Mx4),disp(X'),diary off,...
  117. disp('Press any key to continue.'),pause,diary on,...
  118. clc,disp(Mx5),disp(P),disp(Mx6),disp(L),disp(Mx7),disp((P*B)'),...
  119. disp(Mx8),disp(Y'),diary off,...
  120. disp('Press any key to continue.'),pause,diary on,...
  121. clc,disp(''),disp(Mx9),disp(U),disp(Mx10),disp(Y'),...
  122. disp(Mx11),disp(X'),disp(Mx12),disp(R'),diary off, echo on
  123. pause % Press any key to perform LU factorization with pivoting.
  124.  
  125. clc;
  126. n = 4;                      % Experiment with n=4,5,... keep n small.
  127. A = hilb(n);                % Test the program with a Hilbert matrix.
  128.  
  129. B = ones(n,1);              % Enter B as a column vector.
  130.  
  131. [LU,row,det1] = lufact(A);
  132.  
  133. [X,Y] = lusolv(LU,B,row);
  134.  
  135. R = B - A*X;                % Investigate the residual.
  136.  
  137. pause % Press any key to continue.
  138.  
  139. [n n] = size(A);
  140. I = eye(n);
  141. for k = 1:n,
  142.   P(k,:) = I(row(k),:);
  143. end
  144. L = I;
  145. for c = 1:n-1,
  146.     L(c+1:n,c) = A(row(c+1:n),c);
  147. end
  148. U = zeros(n,n);
  149. for c = 1:n,
  150.     U(1:c,c) = A(row(1:c),c);
  151. end
  152.  
  153. clc,echo off,diary output,...
  154. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(A),disp(Mx3),disp(B'),...
  155. disp(Mx4),disp(X'),diary off,...
  156. disp('Press any key to continue.'),pause,diary on,...
  157. clc,disp(Mx5),disp(P),disp(Mx6),disp(L),disp(Mx7),disp((P*B)'),...
  158. disp(Mx8),disp(Y'),diary off,...
  159. disp('Press any key to continue.'),pause,diary on,...
  160. clc,disp(''),disp(Mx9),disp(U),disp(Mx10),disp(Y'),...
  161. disp(Mx11),disp(X'),disp(Mx12),disp(R'),diary off, echo on
  162.  
  163.  
  164.  
  165.